"""
    VariationalQC

Variational quality control using robust cost functions. Instead of rejecting
outliers, VarQC down-weights them in the cost function using robust M-estimators.

# Mathematical Foundation

## Standard Quadratic Cost
```
J_obs = ½Σᵢ (d_i / σ_i)²
```

## Huber Norm (Robust Cost)
```
ρ(r) = { r²/2           if |r| ≤ c
       { c|r| - c²/2    if |r| > c
```

## Influence Function
```
ψ(r) = dρ/dr = { r      if |r| ≤ c
               { c·sign(r)  if |r| > c
```

## Effective Weight
```
w(r) = ψ(r)/r = { 1        if |r| ≤ c
                { c/|r|    if |r| > c
```

# Implementation

VarQC modifies observation errors adaptively:
```
σ_eff = σ_obs / √w(r)
```

Or equivalently, applies weight in cost function:
```
J = ½Σᵢ w(rᵢ) * rᵢ²
```

# Robust Functions

**Huber**: Linear tail (most common in NWP)
**Bisquare (Tukey)**: Zero weight for large outliers
**Fair**: Softer transition than Huber
**Cauchy**: Very robust but nonconvex
"""
module VariationalQC

using LinearAlgebra
using Statistics
using ...ObservationTypes

# Export types and functions
export VariationalQCConfig, VariationalQCResult
export perform_variational_qc, apply_robust_cost, compute_effective_weight
export huber_norm, bisquare_norm, fair_norm, cauchy_norm

"""
    VariationalQCConfig{T<:AbstractFloat}

Configuration for variational quality control.
"""
struct VariationalQCConfig{T<:AbstractFloat}
    cost_threshold::T
    influence_radius::T
    max_iterations::Int
    convergence_tolerance::T
    robust_function::Symbol  # :huber, :bisquare, :fair, :cauchy
    huber_constant::T
    bisquare_constant::T
    outlier_probability::T

    function VariationalQCConfig{T}(;
        cost_threshold = T(4.0),
        influence_radius = T(200_000),  # 200 km
        max_iterations = 10,
        convergence_tolerance = T(0.01),
        robust_function = :huber,
        huber_constant = T(1.345),      # 95% efficiency for normal
        bisquare_constant = T(4.685),
        outlier_probability = T(0.1)
    ) where T
        new{T}(cost_threshold, influence_radius, max_iterations,
               convergence_tolerance, robust_function, huber_constant,
               bisquare_constant, outlier_probability)
    end
end

VariationalQCConfig() = VariationalQCConfig{Float64}()

"""
    VariationalQCResult{T<:AbstractFloat}

Result of variational QC for a single observation.
"""
struct VariationalQCResult{T<:AbstractFloat}
    passed::Bool
    failure_reason::String
    robust_cost::T
    effective_weight::T
    inflated_error::T
    qc_flag::Int
end

"""
    perform_variational_qc(obs::AbstractObservation{T},
                          background_value::T,
                          config::VariationalQCConfig{T}) where T

Perform variational quality control.
"""
function perform_variational_qc(obs::AbstractObservation{T},
                               background_value::T,
                               config::VariationalQCConfig{T}) where T
    # Compute innovation
    innovation = obs.observed_value - background_value

    # Compute normalized innovation
    total_error = obs.observation_error.total_error
    normalized_innovation = innovation / total_error

    # Apply robust cost function
    robust_cost = apply_robust_cost(normalized_innovation, config)

    # Compute effective weight
    eff_weight = compute_effective_weight(normalized_innovation, config)

    # Compute inflated error (effective weight < 1 means error inflation)
    inflated_error = total_error / sqrt(max(eff_weight, T(0.01)))

    # Decision based on robust cost
    passed = robust_cost <= config.cost_threshold

    qc_flag = if robust_cost <= config.cost_threshold * T(0.5)
        0  # Good
    elseif robust_cost <= config.cost_threshold
        1  # Suspect but acceptable
    else
        2  # Bad - downweighted significantly
    end

    failure_reason = passed ? "" : "high_variational_cost"

    return VariationalQCResult{T}(
        passed, failure_reason, robust_cost, eff_weight,
        inflated_error, qc_flag
    )
end

"""
    apply_robust_cost(r::T, config::VariationalQCConfig{T}) where T

Apply robust cost function to normalized residual.
"""
function apply_robust_cost(r::T, config::VariationalQCConfig{T}) where T
    if config.robust_function == :huber
        return huber_norm(r, config.huber_constant)
    elseif config.robust_function == :bisquare
        return bisquare_norm(r, config.bisquare_constant)
    elseif config.robust_function == :fair
        return fair_norm(r, config.huber_constant)
    elseif config.robust_function == :cauchy
        return cauchy_norm(r, config.huber_constant)
    else
        return r^2 / 2  # Standard quadratic
    end
end

"""
    compute_effective_weight(r::T, config::VariationalQCConfig{T}) where T

Compute effective weight for observation.
"""
function compute_effective_weight(r::T, config::VariationalQCConfig{T}) where T
    abs_r = abs(r)

    if abs_r < T(1e-10)
        return T(1.0)
    end

    if config.robust_function == :huber
        c = config.huber_constant
        return abs_r <= c ? T(1.0) : c / abs_r

    elseif config.robust_function == :bisquare
        c = config.bisquare_constant
        if abs_r <= c
            return (1 - (r/c)^2)^2
        else
            return T(0.0)
        end

    elseif config.robust_function == :fair
        c = config.huber_constant
        return T(1.0) / (1 + abs_r / c)

    elseif config.robust_function == :cauchy
        c = config.huber_constant
        return T(1.0) / (1 + (r/c)^2)

    else
        return T(1.0)  # Standard (no downweighting)
    end
end

# Robust cost functions

"""
    huber_norm(r::T, c::T) where T

Huber norm: quadratic for small residuals, linear for large.
"""
function huber_norm(r::T, c::T) where T
    abs_r = abs(r)
    if abs_r <= c
        return r^2 / 2
    else
        return c * abs_r - c^2 / 2
    end
end

"""
    bisquare_norm(r::T, c::T) where T

Bisquare (Tukey) norm: zero weight beyond threshold.
"""
function bisquare_norm(r::T, c::T) where T
    if abs(r) <= c
        term = 1 - (r/c)^2
        return c^2 / 6 * (1 - term^3)
    else
        return c^2 / 6
    end
end

"""
    fair_norm(r::T, c::T) where T

Fair norm: smooth transition to linear tail.
"""
function fair_norm(r::T, c::T) where T
    abs_r = abs(r)
    return c^2 * (abs_r/c - log(1 + abs_r/c))
end

"""
    cauchy_norm(r::T, c::T) where T

Cauchy norm: very robust but nonconvex.
"""
function cauchy_norm(r::T, c::T) where T
    return c^2 / 2 * log(1 + (r/c)^2)
end

end # module VariationalQC
